function optimization_applications()
% OPTIMIZATION_APPLICATIONS 正定矩阵在优化中的应用
% 
% 展示正定矩阵在各种优化问题中的重要作用

fprintf('=== 正定矩阵在优化中的应用 ===\n\n');

%% 1. 二次规划问题
fprintf('1. 二次规划问题\n');
fprintf('目标：求解标准二次规划问题\n\n');

% 问题：min (1/2)x'Hx + c'x
% 其中 H 是正定矩阵
H1 = [4, 1; 1, 3];
c1 = [-2; -3];

fprintf('二次规划问题: min (1/2)x''Hx + c''x\n');
fprintf('Hessian矩阵 H =\n');
disp(H1);
fprintf('线性项系数 c = [%.1f, %.1f]''\n', c1);

% 检查正定性
eigenvals_H1 = eig(H1);
fprintf('Hessian特征值: [%.4f, %.4f]\n', eigenvals_H1);
fprintf('正定性: %s\n', mat2str(all(eigenvals_H1 > 0)));

if all(eigenvals_H1 > 0)
    % 解析解：x* = -H^(-1)*c
    x1_opt = -H1 \ c1;
    fprintf('最优解: x* = [%.6f, %.6f]''\n', x1_opt);
    
    % 最优值
    opt_value1 = 0.5 * x1_opt' * H1 * x1_opt + c1' * x1_opt;
    fprintf('最优目标函数值: %.6f\n', opt_value1);
    
    % 验证最优性条件：梯度为零
    gradient1 = H1 * x1_opt + c1;
    fprintf('最优点处的梯度: [%.2e, %.2e]''\n', gradient1);
    
    % 可视化目标函数（2D情况）
    if size(H1, 1) == 2
        fprintf('目标函数的等高线分析:\n');
        
        % 计算等高线
        x_range = linspace(x1_opt(1) - 2, x1_opt(1) + 2, 50);
        y_range = linspace(x1_opt(2) - 2, x1_opt(2) + 2, 50);
        [X, Y] = meshgrid(x_range, y_range);
        
        Z = zeros(size(X));
        for i = 1:numel(X)
            x_point = [X(i); Y(i)];
            Z(i) = 0.5 * x_point' * H1 * x_point + c1' * x_point;
        end
        
        fprintf('目标函数值范围: [%.2f, %.2f]\n', min(Z(:)), max(Z(:)));
        fprintf('最优点处的函数值: %.6f\n', opt_value1);
    end
else
    fprintf('Hessian矩阵不是正定的，问题可能无界\n');
end

%% 2. 带等式约束的二次规划
fprintf('\n2. 带等式约束的二次规划\n');
fprintf('目标：求解等式约束二次规划问题\n\n');

% 问题：min (1/2)x'Hx + c'x  s.t. Ax = b
H2 = [2, 1; 1, 2];
c2 = [1; 1];
A2 = [1, 1];  % 约束矩阵
b2 = 1;       % 约束右端

fprintf('约束二次规划问题:\n');
fprintf('min (1/2)x''Hx + c''x  s.t. Ax = b\n');
fprintf('H =\n'); disp(H2);
fprintf('c = [%.1f, %.1f]''\n', c2);
fprintf('A = [%.1f, %.1f]\n', A2);
fprintf('b = %.1f\n', b2);

% 检查正定性
eigenvals_H2 = eig(H2);
fprintf('H的特征值: [%.4f, %.4f]\n', eigenvals_H2);

if all(eigenvals_H2 > 0)
    % KKT系统：[H A'; A 0] * [x; λ] = [-c; b]
    KKT_matrix = [H2, A2'; A2, 0];
    KKT_rhs = [-c2; b2];
    
    fprintf('KKT系统矩阵 =\n');
    disp(KKT_matrix);
    
    % 求解KKT系统
    KKT_solution = KKT_matrix \ KKT_rhs;
    x2_opt = KKT_solution(1:2);
    lambda2_opt = KKT_solution(3);
    
    fprintf('最优解: x* = [%.6f, %.6f]''\n', x2_opt);
    fprintf('拉格朗日乘子: λ* = %.6f\n', lambda2_opt);
    
    % 验证约束满足
    constraint_violation = norm(A2 * x2_opt - b2);
    fprintf('约束违反: %.2e\n', constraint_violation);
    
    % 验证KKT条件
    gradient_lagrangian = H2 * x2_opt + c2 - A2' * lambda2_opt;
    fprintf('拉格朗日函数梯度: [%.2e, %.2e]''\n', gradient_lagrangian);
    
    % 最优值
    opt_value2 = 0.5 * x2_opt' * H2 * x2_opt + c2' * x2_opt;
    fprintf('最优目标函数值: %.6f\n', opt_value2);
else
    fprintf('H不是正定的\n');
end

%% 3. 信赖域方法
fprintf('\n3. 信赖域方法\n');
fprintf('目标：使用正定矩阵求解信赖域子问题\n\n');

% 信赖域子问题：min m(s) = (1/2)s'Bs + g's  s.t. ||s|| ≤ Δ
B3 = [3, 1; 1, 2];  % 正定Hessian近似
g3 = [1; 2];        % 梯度
Delta3 = 1.5;       % 信赖域半径

fprintf('信赖域子问题:\n');
fprintf('min (1/2)s''Bs + g''s  s.t. ||s|| ≤ Δ\n');
fprintf('B =\n'); disp(B3);
fprintf('g = [%.1f, %.1f]''\n', g3);
fprintf('Δ = %.1f\n', Delta3);

% 检查正定性
eigenvals_B3 = eig(B3);
fprintf('B的特征值: [%.4f, %.4f]\n', eigenvals_B3);

if all(eigenvals_B3 > 0)
    % 无约束最优解
    s3_unconstrained = -B3 \ g3;
    s3_unconstrained_norm = norm(s3_unconstrained);
    
    fprintf('无约束最优解: s = [%.4f, %.4f]''\n', s3_unconstrained);
    fprintf('无约束解的范数: %.4f\n', s3_unconstrained_norm);
    
    if s3_unconstrained_norm <= Delta3
        % 无约束解在信赖域内
        s3_opt = s3_unconstrained;
        fprintf('最优解在信赖域内，采用无约束解\n');
    else
        % 需要求解约束问题：(B + σI)s = -g, ||s|| = Δ
        fprintf('无约束解超出信赖域，求解约束问题\n');
        
        % 使用二分法求解σ
        sigma_low = 0;
        sigma_high = 10;
        tol = 1e-8;
        max_iter = 50;
        
        for iter = 1:max_iter
            sigma = (sigma_low + sigma_high) / 2;
            B_sigma = B3 + sigma * eye(2);
            
            try
                s_sigma = -B_sigma \ g3;
                s_sigma_norm = norm(s_sigma);
                
                if abs(s_sigma_norm - Delta3) < tol
                    s3_opt = s_sigma;
                    fprintf('收敛于σ = %.6f，迭代次数: %d\n', sigma, iter);
                    break;
                elseif s_sigma_norm > Delta3
                    sigma_low = sigma;
                else
                    sigma_high = sigma;
                end
                
                if iter == max_iter
                    s3_opt = s_sigma;
                    fprintf('达到最大迭代次数，σ = %.6f\n', sigma);
                end
            catch
                sigma_low = sigma;
            end
        end
    end
    
    fprintf('信赖域最优解: s* = [%.6f, %.6f]''\n', s3_opt);
    fprintf('最优解范数: %.6f (约束: ≤ %.1f)\n', norm(s3_opt), Delta3);
    
    % 目标函数值
    obj_value3 = 0.5 * s3_opt' * B3 * s3_opt + g3' * s3_opt;
    fprintf('目标函数值: %.6f\n', obj_value3);
else
    fprintf('B不是正定的，需要修正\n');
end

%% 4. 牛顿法中的正定性修正
fprintf('\n4. 牛顿法中的正定性修正\n');
fprintf('目标：处理非正定Hessian矩阵\n\n');

% 非正定Hessian矩阵
H4 = [1, 2; 2, 1];  % 不定矩阵
g4 = [1; -1];       % 梯度

fprintf('原始Hessian矩阵 H =\n');
disp(H4);
fprintf('梯度 g = [%.1f, %.1f]''\n', g4);

eigenvals_H4 = eig(H4);
fprintf('H的特征值: [%.4f, %.4f]\n', eigenvals_H4);
fprintf('正定性: %s\n', mat2str(all(eigenvals_H4 > 0)));

if ~all(eigenvals_H4 > 0)
    fprintf('Hessian不是正定的，需要修正\n');
    
    % 修正方法1：添加对角项
    tau1 = abs(min(eigenvals_H4)) + 0.1;
    H4_modified1 = H4 + tau1 * eye(2);
    
    fprintf('方法1：对角修正 H + τI，τ = %.2f\n', tau1);
    fprintf('修正后的特征值: [%.4f, %.4f]\n', eig(H4_modified1));
    
    % 牛顿步
    p4_newton1 = -H4_modified1 \ g4;
    fprintf('修正牛顿步: p = [%.4f, %.4f]''\n', p4_newton1);
    
    % 修正方法2：特征值修正
    [V4, D4] = eig(H4);
    eigenvals_modified = max(diag(D4), 0.1);
    H4_modified2 = V4 * diag(eigenvals_modified) * V4';
    
    fprintf('方法2：特征值修正\n');
    fprintf('修正后的特征值: [%.4f, %.4f]\n', eigenvals_modified);
    
    p4_newton2 = -H4_modified2 \ g4;
    fprintf('修正牛顿步: p = [%.4f, %.4f]''\n', p4_newton2);
    
    % 修正方法3：LDLT分解修正
    fprintf('方法3：LDLT修正\n');
    try
        [L4, D4_ldl, P4] = ldl(H4);
        D4_modified = diag(max(diag(D4_ldl), 0.1));
        H4_modified3 = P4' * L4 * D4_modified * L4' * P4;
        
        fprintf('LDLT修正后的特征值: [%.4f, %.4f]\n', eig(H4_modified3));
        
        p4_newton3 = -H4_modified3 \ g4;
        fprintf('修正牛顿步: p = [%.4f, %.4f]''\n', p4_newton3);
    catch
        fprintf('LDLT分解失败\n');
    end
    
    % 比较不同修正方法
    fprintf('不同修正方法的比较:\n');
    fprintf('方法1步长: %.4f\n', norm(p4_newton1));
    fprintf('方法2步长: %.4f\n', norm(p4_newton2));
    if exist('p4_newton3', 'var')
        fprintf('方法3步长: %.4f\n', norm(p4_newton3));
    end
end

%% 5. 凸优化中的正定性
fprintf('\n5. 凸优化中的正定性\n');
fprintf('目标：分析目标函数的凸性\n\n');

% 二次函数 f(x) = (1/2)x'Hx + c'x + d
H5_convex = [2, 0.5; 0.5, 3];     % 正定矩阵
H5_concave = [-2, 0.5; 0.5, -1];  % 负定矩阵
H5_indefinite = [1, 2; 2, -1];    % 不定矩阵

matrices = {H5_convex, H5_concave, H5_indefinite};
names = {'正定（凸）', '负定（凹）', '不定（非凸非凹）'};

fprintf('不同Hessian矩阵的凸性分析:\n');
for i = 1:length(matrices)
    H = matrices{i};
    eigenvals = eig(H);
    
    fprintf('%s矩阵:\n', names{i});
    fprintf('  特征值: [%.4f, %.4f]\n', eigenvals);
    
    if all(eigenvals > 0)
        convexity = '严格凸';
    elseif all(eigenvals >= 0)
        convexity = '凸';
    elseif all(eigenvals < 0)
        convexity = '严格凹';
    elseif all(eigenvals <= 0)
        convexity = '凹';
    else
        convexity = '非凸非凹';
    end
    
    fprintf('  凸性: %s\n', convexity);
    
    % 对于凸函数，任何局部最小值都是全局最小值
    if all(eigenvals > 0)
        fprintf('  性质: 任何局部最小值都是全局最小值\n');
    elseif all(eigenvals < 0)
        fprintf('  性质: 任何局部最大值都是全局最大值\n');
    else
        fprintf('  性质: 可能存在鞍点\n');
    end
    fprintf('\n');
end

%% 6. 椭球方法
fprintf('\n6. 椭球方法\n');
fprintf('目标：使用椭球方法求解线性规划\n\n');

% 线性规划问题：min c'x  s.t. Ax ≤ b
c6 = [1; 2];
A6 = [1, 1; -1, 1; 0, -1];
b6 = [3; 1; 0];

fprintf('线性规划问题: min c''x  s.t. Ax ≤ b\n');
fprintf('c = [%.1f, %.1f]''\n', c6);
fprintf('A =\n'); disp(A6);
fprintf('b = [%.1f, %.1f, %.1f]''\n', b6);

% 椭球方法的初始椭球
x6_center = [0; 0];  % 椭球中心
E6 = 10 * eye(2);    % 椭球矩阵（正定）

fprintf('初始椭球中心: [%.1f, %.1f]''\n', x6_center);
fprintf('初始椭球矩阵 E =\n');
disp(E6);

% 椭球方法迭代（简化版本）
max_iter_ellipsoid = 10;
tolerance = 1e-6;

fprintf('椭球方法迭代:\n');
fprintf('迭代    中心点              椭球体积      违反约束\n');
fprintf('----    ------              --------      --------\n');

for iter = 1:max_iter_ellipsoid
    % 检查约束违反
    constraint_values = A6 * x6_center - b6;
    violated_constraints = find(constraint_values > tolerance);
    
    ellipsoid_volume = sqrt(det(E6));
    
    fprintf('%4d    [%6.3f, %6.3f]    %8.4f      %d\n', ...
            iter, x6_center, ellipsoid_volume, length(violated_constraints));
    
    if isempty(violated_constraints)
        fprintf('找到可行解\n');
        break;
    end
    
    % 选择最违反的约束
    [~, most_violated_idx] = max(constraint_values);
    a = A6(most_violated_idx, :)';
    
    % 椭球更新公式
    n = length(x6_center);
    
    % 计算切割平面的法向量在椭球度量下的长度
    norm_a_E = sqrt(a' * E6 * a);
    
    if norm_a_E > tolerance
        % 更新椭球中心
        x6_center = x6_center - (1/(n+1)) * (E6 * a) / norm_a_E;
        
        % 更新椭球矩阵
        E6 = (n^2/(n^2-1)) * (E6 - (2/(n+1)) * (E6 * a * a' * E6) / (a' * E6 * a));
        
        % 确保正定性
        eigenvals_E6 = eig(E6);
        if any(eigenvals_E6 <= 0)
            fprintf('椭球矩阵失去正定性，算法终止\n');
            break;
        end
    else
        fprintf('椭球退化，算法终止\n');
        break;
    end
end

% 最终结果
final_obj_value = c6' * x6_center;
fprintf('椭球方法最终解: x = [%.4f, %.4f]''\n', x6_center);
fprintf('目标函数值: %.4f\n', final_obj_value);

%% 7. 半正定规划
fprintf('\n7. 半正定规划\n');
fprintf('目标：展示半正定约束的优化问题\n\n');

% SDP问题：min <C, X>  s.t. <A_i, X> = b_i, X ⪰ 0
% 简单例子：2×2矩阵变量
fprintf('半正定规划问题:\n');
fprintf('min <C, X>  s.t. <A_i, X> = b_i, X ⪰ 0\n');

% 目标矩阵
C7 = [1, 0; 0, 1];
fprintf('目标矩阵 C =\n');
disp(C7);

% 约束矩阵
A7_1 = [1, 0; 0, 0];  % 约束1：X(1,1) = 1
A7_2 = [0, 1; 1, 0];  % 约束2：X(1,2) + X(2,1) = 0
b7 = [1; 0];

fprintf('约束矩阵 A1 =\n'); disp(A7_1);
fprintf('约束矩阵 A2 =\n'); disp(A7_2);
fprintf('约束右端 b = [%.1f, %.1f]''\n', b7);

% 可行解示例
X7_feasible = [1, 0; 0, 2];  % 满足约束的正定矩阵
fprintf('可行解示例 X =\n');
disp(X7_feasible);

% 验证约束
constraint1_value = trace(A7_1' * X7_feasible);
constraint2_value = trace(A7_2' * X7_feasible);
fprintf('约束验证: <A1, X> = %.1f (应为1)\n', constraint1_value);
fprintf('约束验证: <A2, X> = %.1f (应为0)\n', constraint2_value);

% 验证正定性
eigenvals_X7 = eig(X7_feasible);
fprintf('X的特征值: [%.4f, %.4f]\n', eigenvals_X7);
fprintf('正定性: %s\n', mat2str(all(eigenvals_X7 > 0)));

% 目标函数值
obj_value7 = trace(C7' * X7_feasible);
fprintf('目标函数值: %.4f\n', obj_value7);

%% 总结
fprintf('\n=== 优化应用总结 ===\n');
fprintf('1. 二次规划：正定Hessian保证凸性和唯一解\n');
fprintf('2. 约束优化：KKT系统的求解依赖正定性\n');
fprintf('3. 信赖域方法：正定模型确保下降方向\n');
fprintf('4. 牛顿法：非正定Hessian需要修正\n');
fprintf('5. 凸优化：正定性决定函数的凸性\n');
fprintf('6. 椭球方法：椭球矩阵必须保持正定\n');
fprintf('7. 半正定规划：正定性作为约束条件\n\n');

fprintf('正定矩阵在优化中的重要性：\n');
fprintf('• 保证问题的凸性\n');
fprintf('• 确保算法的收敛性\n');
fprintf('• 提供数值稳定性\n');
fprintf('• 简化理论分析\n\n');

fprintf('实际应用中的注意事项：\n');
fprintf('• 检查矩阵的正定性\n');
fprintf('• 必要时进行正定性修正\n');
fprintf('• 考虑数值精度问题\n');
fprintf('• 选择合适的算法\n\n');

end